{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pylab as pl\n",
    "\n",
    "pl.rcParams[\"figure.figsize\"] = 4, 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gmm_base import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lamb = 0.5\n",
    "xmax = 4 / lamb\n",
    "sample_size = int(1e4)\n",
    "num_clusters = 3\n",
    "epochs = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.RandomState(0).exponential(1.0 / lamb, sample_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Gaussian_pdf(ex):\n",
    "    return np.exp(\n",
    "        -0.5 * (ex - x.mean()) ** 2 / x.std() ** 2\n",
    "        - 0.5 * np.log(2 * np.pi * x.std() ** 2)\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.4191023651728891, 8.0)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAD4CAYAAADhGCPfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVxVdf7H8ddHFlHcSnFDE3TMRJFUxAXFXSkdLa0mbZHKscatbMpsmaaxZdRpoqacmjKzxcatdJyyzKX0p5mgprllmmmiTqIGqYmyfH9/HEDAA1zgwrkXPs/Hg4fcc88993MV35zzPd9FjDEopVRB1ZwuQCnlmTQclFK2NByUUrY0HJRStjQclFK2fJ164wYNGpiQkBCn3l6pKmvr1q0njTFBxe3nWDiEhISwZcsWp95eqSpLRA67sp9eViilbGk4KKVsaTgopWw51uagPF96ejpJSUmkpaU5XYoqhYCAAJo1a4afn1+pXq/hoAqVlJRE7dq1CQkJQUScLkeVgDGGU6dOkZSURGhoaKmOoZcVqlBpaWnUr19fg8ELiQj169cv01mfhoMqkgaD9yrrv52Gg1LKloaD8mg//fQTo0ePpmXLlnTu3Jnu3buzdOnScn3PLVu2MHny5HJ9D2/gFQ2S0TPWcjTlfL5twfVqsHFaP4cqUhXBGMMNN9zAmDFjeP/99wE4fPgwy5cvL9f3jYyMJDIyslzfwxt4xZnD0ZTzHJoxJN9XwbBQlc/atWvx9/fnvvvuy93WokULJk2axKFDh+jVqxedOnWiU6dOfPnllwB88cUXDB06NHf/iRMnMm/ePACmTZtGWFgYHTp04KGHHgJg8eLFtG/fnoiICGJiYi47RkJCAj169KBjx4706NGDffv2ATBv3jxGjBhBbGwsrVu3ZurUqeX+91HRXDpzEJFY4CXAB5hjjJlR4Pk44G/A0exNrxhj5rixTuUJ3hpy+bZ2N0DU7+HirzD/5sufv3Y0dLwNzp2CRXfmf+6uj4t8u927d9OpUyfb5xo2bMiqVasICAhg//79jBo1qsixOqdPn2bp0qV8++23iAgpKSkATJ8+nZUrVxIcHJy7La9rrrmG9evX4+vry+rVq3nsscf44IMPANi+fTtff/011atXp02bNkyaNInmzZsX+Zm8SbHhICI+wGxgIJAEJIrIcmPMngK7LjTGTCyHGpUCYMKECWzYsAF/f39Wr17NxIkT2b59Oz4+Pnz33XdFvrZOnToEBAQwduxYhgwZkntmEB0dTVxcHLfccgsjRoy47HWpqamMGTOG/fv3IyKkp6fnPte/f3/q1q0LQFhYGIcPH65a4QBEAQeMMQcBRGQBMBwoGA6qsivqN71/zaKfD6xf7JlCQe3atcv9LQ0we/ZsTp48SWRkJPHx8TRq1IgdO3aQlZVFQEAAAL6+vmRlZeW+Juc+v6+vLwkJCaxZs4YFCxbwyiuvsHbtWl577TU2b97Mxx9/zLXXXsv27dvz1fCnP/2Jvn37snTpUg4dOkSfPn1yn6tevXru9z4+PmRkZJTo83k6V9ocgoEjeR4nZW8raKSIfCMiS0TENj5FZJyIbBGRLcnJyaUoV1Ul/fr1Iy0tjVdffTV326+//gpYv9GbNGlCtWrVePfdd8nMzASsNok9e/Zw4cIFUlNTWbNmDQBnz54lNTWV66+/nhdffDE3BL7//nu6du3K9OnTadCgAUeOHMlXQ2pqKsHB1o97TttFVeFKONj1pCg4n/1/gRBjTAdgNfC23YGMMa8bYyKNMZFBQcXONaGqOBFh2bJlrFu3jtDQUKKiohgzZgwzZ85k/PjxvP3223Tr1o3vvvuOwMBAAJo3b84tt9xChw4duO222+jYsSMAZ86cYejQoXTo0IHevXsTHx8PwMMPP0x4eDjt27cnJiaGiIiIfDVMnTqVRx99lOjo6NwAqjKMMUV+Ad2BlXkePwo8WsT+PkBqccft3LmzcVWLRz5yaZtyrz179jhdgioju39DYIsp5v+nMcalM4dEoLWIhIqIP3ArkO9Gs4g0yfNwGLC3rKGllHJWsQ2SxpgMEZkIrMQ6K5hrjNktItOxEmg5MFlEhgEZwGkgrhxrVkpVAJf6ORhjVgArCmx7Ms/3j2JdbiilKgmv6CGplKp4Gg5KKVsaDkopW14xKlN5BrvRsWXh1MjaHj165A7UKslzRXnqqaeoVatW7oCuykDDQbksZ3Ssu4RMK1l3anex+8+fmZmJj49PqYKhstLLCuXR3nvvPaKiorj22mu59957yczMpFatWjzyyCN07tyZAQMGkJCQQJ8+fWjZsmXuXA/z5s1j+PDhxMbG0qZNG/7yl7/kHrNWrVqANTS7b9++jB49mvDw8HzPAcyaNYvw8HAiIiKYNm0aAG+88QZdunQhIiKCkSNH5nbnrow0HJTH2rt3LwsXLmTjxo25oy/nz5/PuXPn6NOnD1u3bqV27do88cQTrFq1iqVLl/Lkk7l32ElISGD+/Pls376dxYsX2w7pTkhI4Nlnn2XPnvzjCD/55BOWLVvG5s2b2bFjR+58DSNGjCAxMZEdO3bQtm1b3nzzzfL9S3CQXlYoj7VmzRq2bt1Kly5dADh//jwNGzbE39+f2NhYAMLDw6levTp+fn6Eh4dz6NCh3NcPHDiQ+vXrA9Z/6g0bNlw2w1NUVJTt1O2rV6/mrrvuombNmgBceeWVAOzatYsnnniClJQUzp49y+DBg93+uT2FhoPyWMYYxowZw1//+td8259//vncmZWrVauWO3S6WrVq+YZNF5x92W425pwBW3bvbbd/XFwcy5YtIyIignnz5vHFF1+U6DN5E6+9rAiuV4OQaR/nfkXPWOt0ScrN+vfvz5IlSzhx4gRgzeZ0+LBLC0QDsGrVKk6fPs358+dZtmwZ0dHRLr920KBBzJ07N7dN4fTp04A1urNJkyakp6czf/78Enwa7+O1Zw4Fb4E51fJdleQEsjuPV5SwsDCeeeYZBg0aRFZWFn5+fsyePdvl4/fs2ZM77riDAwcOMHr06BJNGhsbG8v27duJjIzE39+f66+/nueee46nn36arl270qJFC8LDwzlz5ozLx/Q6rgzdLI+vsg7ZLs0+qmS8ecj2W2+9ZSZMmOB0GY4r7yHbSqkqyGsvK5QqSlxcHHFxcU6X4dX0zEEpZUvDQSllS8NBKWVLw0EpZUvDQXk0J1bZBs9daTsuLo4lS5ZUyHtpOCi3ODVnDue+2pxv27mvNnNqTumXTDXZq2zHxMRw8OBBtm7dyoIFC0hKSiprucWKjIzkH//4R7m/jyfTcFBlkvzyKwAEtA/n6JQpuQFx7qvNHJ0yhYD24fn2K4miVtkGvGal7ZCQEB577DG6d+9OZGQk27ZtY/DgwbRq1YrXXnsNsILw4Ycfpn379oSHh7Nw4cLc7RMnTiQsLIwhQ4bkdiUH2Lp1K71796Zz584MHjyY48ePl/jvuCjaz0GVycnZs/k1IQEA34YN+XHsWHwbNiTjxAmqt2rFydmzrX0SEwmaVLJ1lotaZRu8a6Xt5s2bs2nTJqZMmUJcXBwbN24kLS2Ndu3acd999/Hhhx+yfft2duzYwcmTJ+nSpQsxMTFs2rSJffv2sXPnTn766SfCwsK4++67SU9PZ9KkSfznP/8hKCiIhQsX8vjjjzN37twS/R0XRcNBuY1PnTpWMBw7hm/TpvjUqePW4+ddZTsxMZH09HSvWWl72LBhgDXE/OzZs9SuXZvatWsTEBBASkoKGzZsYNSoUfj4+NCoUSN69+5NYmIi69evz93etGlT+vWzxhTt27ePXbt2MXDgQMCayapJkyaXvW9ZaDioMmvx7jvApUuJBuP/wM//XkCDCRMI7NYVgL3XtC3xcYtaZRvwqpW28w4rz/uanGHm1pAHe3ZDx40xtGvXjk2bNhX6urLSNgdVJg0mTAAuBUNwfDxBkycTHB+frw0iZ7+SKGqVbahcK23HxMSwcOFCMjMzSU5OZv369URFRRETE8OCBQvIzMzk+PHjfP755wC0adOG5OTk3HBIT09n9+7dbq1JzxxUmeS0I6Tt2klwfHzumUJgt64Ex8eTtmsngd26lri9AS6tsj1lyhRmzZpFUFAQgYGBzJw5E4Dx48czcuRIFi9eTN++fW1X2m7dunW+lbaHDx9OWloaxph8K23v378fYwz9+/cnIiKCdevW5dYxdepUxowZwwsvvJB7Wu9uN954I5s2bSIiIgIRYdasWTRu3Jgbb7yRtWvXEh4eztVXX03v3r0B8Pf3Z8mSJUyePJnU1FQyMjJ44IEHaNeundtqkqJOZ8pTZGSkKarxKK+QaR8XO+uxK/uoktm7dy9t25b8ckB5Drt/QxHZaowpdnILvaxQStnScFBK2dJwUEVy6rJTlV1Z/+00HFShAgICOHXqlAaEFzLGcOrUqdzbu6Xh0t0KEYkFXgJ8gDnGmBmF7HcTsBjoYoxxrbVReaxmzZqRlJREcnKy06WoUggICKBZs2alfn2x4SAiPsBsYCCQBCSKyHJjzJ4C+9UGJgObLz+K8kZ+fn62C76oqsGVM4co4IAx5iCAiCwAhgN7Cuz3NDALcGSZ4YLTpju1grNSlYUr4RAM5O0ylgR0zbuDiHQEmhtjPhIRR8JB17FQyr1caZC8vGM35LZQiUg1IB74Y7EHEhknIltEZItexyrl2VwJhyQg7zCzZsCxPI9rA+2BL0TkENANWC4il/XAMsa8boyJNMZEBgUFlb5qpVS5cyUcEoHWIhIqIv7ArcDynCeNManGmAbGmBBjTAjwFTBM71Yo5d2KDQdjTAYwEVgJ7AUWGWN2i8h0ERlW3gUqpZzhUj8HY8wKYEWBbU8Wsm+fspellHJapR2yrbc2lSqbShsOemtTqbLRsRVKKVsaDkopWxoOSilbGg5KKVsaDkopWxoOSilbGg5KKVsaDkopWxoOSilbGg5KKVsaDkopWxoOSilbGg5KKVuVdlRmQQWHcOds02HcStmrMuFgFwI6jFupwullhVLKloaDUsqWhoNSypaGg1LKloaDUsqWhoNSypaGg1LKloaDUsqWhoNSypaGg1LKVpXpPu2K6BlrOZpyPvexjr1QVZmGQx5HU85zaMaQ3Mc69kJVZXpZoZSypeGglLJVpS8rCs7xEFyvhoPVKOVZXAoHEYkFXgJ8gDnGmBkFnr8PmABkAmeBccaYPW6u1e20sVGpwhV7WSEiPsBs4DogDBglImEFdnvfGBNujLkWmAW84PZKlVIVypU2hyjggDHmoDHmIrAAGJ53B2PML3keBgLGfSUqpZzgymVFMHAkz+MkoGvBnURkAvAg4A/Ynq+LyDhgHMBVV11V0lqVUhXIlTMHsdl22ZmBMWa2MaYV8AjwhN2BjDGvG2MijTGRQUFBJatUKVWhXAmHJKB5nsfNgGNF7L8AuKEsRSmlnOdKOCQCrUUkVET8gVuB5Xl3EJHWeR4OAfa7r0SllBOKbXMwxmSIyERgJdatzLnGmN0iMh3YYoxZDkwUkQFAOvAzMKY8i1ZKlT+X+jkYY1YAKwpsezLP9/e7uS6llMO0+7RSypaGg1LKVpUeW1EcXV9TVWUaDkXQ9TVVVaaXFUopWxoOSilbGg5KKVsaDkopWxoOSilbGg5KKVsaDkopW9rPoYTsJqXVTlGqMtJwKKGCQaCdolRlpZcVSilbGg5KKVsaDkopWxoOSilbGg5KKVsaDkopWxoOSilbGg5KKVsaDkopWxoOSilb2n3azaJnrOVoyvl823T8hfJGGg5lZDcQ69CMIfn20fEXyhtpOJSRnhGoykrbHJRStjQclFK2NByUUrY0HJRStjQclFK2XAoHEYkVkX0ickBEptk8/6CI7BGRb0RkjYi0cH+p3u/UnDmc+2pzvm3nvtrMqTlzHKpIqcIVGw4i4gPMBq4DwoBRIhJWYLevgUhjTAdgCTDLnUXetnelOw/nmID24RydMiU3IM59tZmjU6YQ0D7c4crcI/nlV5wuQbmRGGOK3kGkO/CUMWZw9uNHAYwxfy1k/47AK8aY6KKOGxkZabZs2eJSkXuvaUvNLl1c2tfzZLH96P/oEhqIj89FzifD+f1J+DWsT0ZyMoHXNMKnTh0y0gMxxsfpYsvk18RE2n671+kyVDFEZKsxJrK4/VzpBBUMHMnzOAnoWsT+9wCfFFLUOGAcwFVXXeXCW3sbg0gGxvjh4/srQY0S8K9+hhatsnL3OOUXTnpqQ9KPHaNBuzMEtUuyXmkgIz2Q0yfDSTvfCMhExGCM9lNTznDlJ09sttmebojI7UAk0NvueWPM68DrYJ05uFgjAC3efacku1eczHQ4vBH2/hf2fgThN8HgZyE9Df59KzTpwN82p/H9rzX42dTGP+0ik04tJXTcPfy8cBE1b7uLwNZByMn9+CXvpVHPKdC0I+xeCkvvg5Be0HogtLkO6nl2oO69pq3TJSg3ciUckoDmeR43A44V3ElEBgCPA72NMRfcU57lvTYDedadB3SXlY/D1+9BWgr41YTfDIBWfa3n/ALgzmUAPDzQ2pTTxvBQ5O0sf/B+avboxdEpUwiOjyew7035j12/NXS+Cw6sgk+mwiePWMe+6S2oUa8CP6TrGkyY4HQJyo1cCYdEoLWIhAJHgVuB0Xl3yG5n+BcQa4w54e4i57cd7BnhkJUFRzZDi+7W44tnrd/qYTdAq37gX7PIl6ft2klwfDzfLDsJQGC3rgTHx5O2ayeB3QpcqTVuD9fNAGbAqe9h52LrvQPqWs/v+xSaRUJgAzd/yNILmjTR6RKUGxUbDsaYDBGZCKwEfIC5xpjdIjId2GKMWQ78DagFLBYRgB+NMcPKse6KZQzs+Q98/iyc/A7uXQ9NImDoiyB2V1326o8da32z7NIozcBuXS8Phste2Ar65LmDnHERlt5rXdJE/R56TIbA+iX5REoVy6XWLmPMCmBFgW1P5vl+gJvr8hwH18Hqp+DYNghqCyPfhIbZd3JLEAxu5esP96yCdTNh40uQOAe63muFhIdecijvoz0ki3LhDCy6A86egOH/hD9stBocffycrgyCroab3oTxX1mXNhvi4dQBp6tSlYjeJysoPQ12vA+d4qB6bbhjKTRsZzUweqKG18DN8yA1Ceo2s7Zt+ieExljtFkqVkp455HVwHbzaAz6aAofWW9uCO3tuMOSVEwxpqbDhBfhXDHz6GKSfL/p1ShVCwwHgwllYPgneGQYmyzpbaNnH6apKJ6AuTEiATnfCV7OtkDi6zemqlBfScABYMBq2vQvR98P4TdZtSW9W80r47YtWyF04C+8Mt84olCqBqtvmkJVpnSX4+EGfR6H3VAjp6XRV7tWqH4z/Eo59fal/xLmTHtU3QnmuqnnmcOZ/8PZvYe0z1uMW3StfMOSoccWlM6Ht/4Z/dLK6eitVjKp35nD4S1gcZ92m7HSnIyXYTWdfIbNYt+gB9VvCwtuh+0QY8JRn3JZVHqnqhIMx8NU/4bM/wZWhcMcyaFRwWoqKUTAIKmxdiytawN0r4bMnYNMrcCTBug1aN7hi3l95lapzWfHzD7D6L9boxt9/7lgwOM63Olz/N2sA14m9cNS1OTVU1VP5zxzSfoGAOnBlSxj3udX12aluz56k/QgI7X1pTMap760xHEplq9zhcCTRuk056GmIuBUatXO6IpcVXHOzXNolcoLh+A54oz9EjYOB08Gncv9YKNdU3p+CbxbBfyZCnSbW5CkerGADZc62vGtulmu7RMMwiLzb6jSV/K3VDhFQp/zeT3mFyhcOxsDnz8H6WdCiJ/zuXatTkAdzfL1NHz+4fpbVDvPRgzA3Fm5bdKlLtqqSKl+D5OGNVjB0vN3qIejhweBROsfBbYsh5Uf4ZqHT1SiHVZ4zB2OshsaQnhD3MbSI1obH0vhNf/jDBqibPV/lxXPgH+hsTcoRlePMIeUIvNHPum8PVkBoMJTeFSFQrZp1BvFyZ9g6z+mKlAO8Pxx+2g1vDrQmOtHhye5V40po1B7+ez98MdM6O1NVhneHw6GNMPc6QODuT6Gl7Yz4qrSq14JR/4aI0fDFc/DRA5CZ4XRVqoJ4b5vDse3w7o1Wl+DbP4R6zYt/jSo5Hz+44Z9Qu7E1iUzd5hDzkNNVqQrgveHQONyaf6HbH6rEHQm7vhB2+5TLbVERGPBnqxNZm+vcf3zlkbwrHIyBr16FsOHWYKF+jztdUYVx5T99uQ/gCs9eeOfCGVgx1QqM2o3L9z2VY7ynzSEry1r5aeWjsO1tp6up2pL3Wet4vDnQGpOhKiWvCAc/MuDDsZDwujUPQe9pxb9IlZ9mkRD3X6sPxJuDrPYfVel4/mXFxXPM8Xsedn0DA/4CPR9wuiKvUa6Dt4I7W3NDvHsjzBtq9azMWSZQVQqeHw6Z6dSXX2DYK9DpDqer8SpHU86X7+CtBq3hns9g2XhrgJuqVDz/sqJGPYZffFqDwVPVaWqtJn5FiNUu9MP/OV2RchPPDwcgEx+nS1Cu2PY2vD0UNrzodCXKDbwiHJSXuPY2aDcCVv/ZmqcyK8vpilQZeH6bg3KZ3azWFcrX31qFPLABfPmytUbGsJd1hmsvpeFQiTg+aQxYozmvmwWBQbD+eeh6HzS91umqVCm4dFkhIrEisk9EDojIZZ0MRCRGRLaJSIaI3OT+MpVXEbFWEJuYeCkYMi46W5MqsWLPHETEB5gNDASSgEQRWW6M2ZNntx+BOEBH5HiwwuaqLLczjitaWH/uWAgb4q2+EDpAzmu4clkRBRwwxhwEEJEFwHAgNxyMMYeyn9MWKA9mFwIVsqBOnabwyzGrN+XtS7xqFvCqzJXLimDgSJ7HSdnbSkxExonIFhHZkpycXJpDKG8U2gvu/gQw1vwb2hfCK7hy5mA331qppgQyxrwOvA4QGRmp0wp5gApbt7NRO7hnFbw3Et4bARMSrGUJlcdyJRySgLwXis2AY+VTjqpoFbpuZ73m1oxd336kweAFXLmsSARai0ioiPgDtwLLy7csVWnVvPLS6uZHt8In03TqOQ9VbDgYYzKAicBKYC+wyBizW0Smi8gwABHpIiJJwM3Av0Rkd3kWrSqJg+tg86uwYBRcOOt0NaoAlzpBGWNWACsKbHsyz/eJWJcbSrmu14NQ4wr4+I/wViyMWmjN8KU8go6tUM6KvAtGL4LTh+CNvnD6B6crUtm0+7QqUsEJY6Ac7mi0HgBjV8Hmf0G9q9x3XFUmGg6qSAUnjIFyuqPRsC38Nnuo9y/HYcf7ED3FGquhHKF/88rz7FwEa6ZbDZXnU5yupsrScFCep8dkuP55OLAaXu9jLXmoKpxeVqh8HJ8TAqxRnVG/h8YdYNGdMGeA1WgZ2qvia6nCNBxUPqVpaHRllutSzYR9VVe4d701s1STDiWuS5WNhoMqM1dmuS71TNi1G8GNr1nfp6fBxw9aa3Ve2bJMNaviaZuD8h4n91njMl6LgZ1LnK6m0tMzB1VijrVLNImA+zbAknvgg3vgh/UQ+1fwD6yY969iNBxUiTk6V2W9q+CuFfD5s9YU+GkpcMs7ztVTiWk4KO/j4wcDnoLfDICa9a1tF86AT3VrBmzlFtrmoLxXSE+rZyXARw/CG/3g+DfO1lSJaDioyqH9CDj7k9VpatWTcPFXpyvyenpZodyusFmuy1Wb62BighUMG1+C3cvglrehacfyfd9KTMNBuZ0rDZZ2AWK3T95jFduRqsYV1gpbHX5nhUStxtZ2Y6xel6pENByUI1wJkOgZay+7ZepSR6qQnjB2jRUIxlgT2l7VA3pMBD8HuoN7KQ0H5bGKC5AiF+nJOVO4eA6q14bPn4Ft78Cg6RB2g55JuEDDQXktlxbpqV7L6gfxw//Bp9NgcZzVDnHzPLgipCLK9FoaDqpSKXQdjtBeMG4dfLMQvn73UnvE+Z+ttgp1GQ0HVakUuQ6Hjy/RK5twNGUS/GkN/qSzxv8h/ufblC63T4fQ3nq5kYeGg6pS8o0OTU+DhMNU/+wFeGe4dbnRcwpcMxSq+ThbqAfQTlCq6vILgOj76XXhRRj6ojUl3aI74Yd1TlfmEfTMQVV5DerVJWSJP9V4mgHVtrJnURYbHgW+mAGpSdBlLDS91ukyK5yGg6ry8rdT/PZSO0VGGuz6wGrAbNQeIm6F8JuhdmNH6qxoGg6qUivT3BMDnoLoB3j+78/S//haOv70BB988il/TB9PcN0ANj7QqUR3Oko1VZ6DNBxUpVbm/3w16vHK2T48NONvcHI/IxFGNvgNdz43l4wZN7EpK4xPs6JYldmZExQdFC738PQQGg5KFVDo2UaD1rnb3hk/ELbcT689y+l1ei7P+s2Fhu3g5rcgqE1Fl1wuNByUKsCls416za3Ljv5/hhN7Yf9ncPBzqJO9EPCGePh+LbSIhhY9IDgS/GvmO0ShHbY8hIaDUmUhAo3CrK+eD1za7l/LujX6xQzAQDVfuKo7xH1kPf/LMTY+3Mua1Sqbp11maDgoVR6ifm99nU+BIwnw45eQmX7p+fdvgeR90OBqaGiFS6RkAfnXJbVbyDiv8jzb0HBQqjzVqAdXD7K+8ur5IBzfDj/tgcMbYeci7q3Rk5Bp1wCwzP8Jfja1edC/CSN/2wPqNrdm367fKt9h7M423HVXxKVwEJFY4CXAB5hjjJlR4PnqwDtAZ+AU8DtjzKESV6NUVdF+hPWV43wKAy+e41DdYMi4CEsXw6kDkLIJVn1q7dNjMgx62ppM9+VIqN2Y92r6suiJf3GKOqzO7MRW04bQutU4NL6BFUw1ruDqZ74sVYnFhoOI+ACzgYFAEpAoIsuNMXvy7HYP8LMx5jciciswE/hdqSpSqiqqUc/6AmsG7ZvnXXouLRVSjljzUoB1edJ6AJz5iZ4N/wfnjsG5ZP4wqCd0GwInvoV/ds19+XcBcOHPvjyWPpYPsmJcLsmVM4co4IAx5iCAiCwAhgN5w2E48FT290uAV0REjDHG5UqUUvYC6kLjupce17wShs/Ov48xkJVpfV83GG7/wGrvOP8zXPiF6mm/8Pew4fw9uBMy07W3dSUcgoEjeR4nAV0L28cYkyEiqUB94GTenURkHDAu++FZEdnnWpk0kJn5j1UJNAD9TLAg+PkAAAL4SURBVF6gEn2m6TnfuNQRw5VwsBvgXvCMwJV9MMa8DrzuwnvmP7jIFmNMZElf58n0M3mHyvqZXNnPlSHbSUDzPI+bAccK20dEfIG6wGlXClBKeSZXwiERaC0ioSLiD9wKLC+wz3JgTPb3NwFrtb1BKe9W7GVFdhvCRGAl1q3MucaY3SIyHdhijFkOvAm8KyIHsM4YbnVznSW+FPEC+pm8Q5X9TKK/4JVSdnSaOKWULQ0HpZQtjw8HEYkVkX0ickBEpjldT1mJSHMR+VxE9orIbhG53+ma3EVEfETkaxH5yOla3EFE6onIEhH5Nvvfq7vTNZWViEzJ/rnbJSL/FpGAwvb16HDI03X7OiAMGCUiYc5WVWYZwB+NMW2BbsCESvCZctwP7HW6CDd6CfjUGHMNEIGXfzYRCQYmA5HGmPZYNxgKvXng0eFAnq7bxpiLQE7Xba9ljDlujNmW/f0ZrB+4YGerKjsRaYY13niO07W4g4jUAWKw7sRhjLlojElxtiq38AVqZPdHqsnlfZZyeXo42HXd9vr/SDlEJAToCGx2thK3eBGYCmQ5XYibtASSgbeyL5XmiEig00WVhTHmKPA88CNwHEg1xnxW2P6eHg4udcv2RiJSC/gAeMAY84vT9ZSFiAwFThhjtjpdixv5Ap2AV40xHYFzgFe3eYnIFVhn3qFAUyBQRG4vbH9PDwdXum57HRHxwwqG+caYD52uxw2igWEicgjr0q+fiLznbElllgQkGWNyzuqWYIWFNxsA/GCMSTbGpAMfAj0K29nTw8GVrtteRUQE6zp2rzHmBafrcQdjzKPGmGbGmBCsf6O1xphCfyN5A2PM/4AjIpIzgrE/+acp8EY/At1EpGb2z2F/imhk9ehp4grruu1wWWUVDdwB7BSR7dnbHjPGrHCwJmVvEjA/+xfTQeAuh+spE2PMZhFZAmzDumv2NUV0pdbu00opW55+WaGUcoiGg1LKloaDUsqWhoNSypaGg1LKloaDUsqWhoNSytb/Aw1yc53lgE8OAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ex = np.linspace(0, xmax)\n",
    "pl.hist(x, ex, density=True, fill=False, histtype=\"step\", label=\"empirical\")\n",
    "\n",
    "pl.plot(ex, Gaussian_pdf(ex), \"--\", label=\"Gaussian\")\n",
    "pl.errorbar(\n",
    "    [x.mean()],\n",
    "    Gaussian_pdf([x.mean()]),\n",
    "    None,\n",
    "    [x.std()],\n",
    "    \"x\",\n",
    "    color=\"C3\",\n",
    "    capsize=2,\n",
    "    label=\"Gaussian mode\",\n",
    ")\n",
    "\n",
    "pl.legend(loc=\"upper right\")\n",
    "pl.xlim(xmax=xmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loglik=-2.111 elapsed=0.0s\n",
      "2 loglik=-2.089 elapsed=0.0s\n",
      "4 loglik=-1.952 elapsed=0.0s\n",
      "8 loglik=-1.887 elapsed=0.1s\n",
      "16 loglik=-1.799 elapsed=0.1s\n",
      "32 loglik=-1.778 elapsed=0.3s\n",
      "64 loglik=-1.775 elapsed=0.8s\n",
      "100 loglik=-1.774 elapsed=1.1s\n"
     ]
    }
   ],
   "source": [
    "model = GMMModel(x[:, None], num_clusters=num_clusters)\n",
    "trainer = GMMTrainer(model)\n",
    "\n",
    "for t, epoch in elapsed(range(100)):\n",
    "    trainer(x[:, None])\n",
    "    if (\n",
    "        np.allclose(np.log2(epoch + 1), np.round(np.log2(epoch + 1)))\n",
    "        or epoch + 1 == 100\n",
    "    ):\n",
    "        loglik = model(mx.nd.array(x[:, None]))[0].mean().asscalar()\n",
    "        print(f\"{epoch+1} loglik={loglik:.3f} elapsed={t:.1f}s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inferred lamb=0.425\n"
     ]
    }
   ],
   "source": [
    "hat_lambda = infer_lambda(model, xmin=0, xmax=xmax)\n",
    "print(f\"inferred lamb={hat_lambda:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mixture_pdf(ex):\n",
    "    log_marg = model(mx.nd.array(ex, dtype=\"float32\"))[0]\n",
    "    return log_marg.exp().asnumpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_prior_ = model.log_prior_.data().asnumpy()\n",
    "mu_ = model.mu_.data().asnumpy()\n",
    "kR_ = model.kR_.data().asnumpy()\n",
    "cov_ = np.linalg.inv(kR_.swapaxes(1, 2) @ kR_)\n",
    "s2_ = np.array([np.diag(c) for c in cov_])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.4, 8.0)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQcAAAD4CAYAAADhGCPfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVhV1frA8e/iMAmKKDgCCikqCoiIU86aQ2pYWlmhOaSZqU23wW7da3Naem2wbpqZ5ZClTeZUmvVzuE445DwiKZiJoCLz9P7+OIAMBzjgwcPB9Xme88jee+193iPwsoe13qVEBE3TtKLsrB2ApmlVk04OmqaZpJODpmkm6eSgaZpJOjlommaSvbXe2NPTU3x9fa319pp2y9qzZ88lEalXVjurJQdfX18iIyOt9faadstSSv1pTjt9WaFpmkk6OWiaZpJODpqmmWS1ew6adWVmZhITE0NaWpq1Q9EqibOzM97e3jg4OFRof50cblExMTHUqlULX19flFLWDkezMBEhPj6emJgY/Pz8KnQMfVlxi0pLS8PDw0MnhmpKKYWHh8cNnRnq5HAL04mhervR769ODpqmmaSTg1alrVq1ihkzZpTaZtGiRZw/f/4mRWTa77//zpAhQwBIT0/njjvuICQkhK+//tqqcd0Im7gh2XXGJmKvpOYve7nXYNu0PlaMSLtZwsPDCQ8PL7XNokWLCAwMpHHjxmYfNysrC3v7yvnx37dvH5mZmezfv79Sjn+z2MSZQ+yVVKJnDM5/FUwUmm2Kjo6mVatWjB8/nsDAQCIiIti4cSNdu3bF39+fXbt2AcZf/ClTpgAwdOhQvvzySwDmzZtHREQEK1euJDIykoiICEJCQkhNTcXX15dLly4BEBkZSa9evQB45ZVXePTRR+nfvz8PP/ww2dnZPPfcc3To0IHg4GDmzZtXYpyjR48mODiYe++9l5SUFADWr19Pq1at6NatG9999x0AFy9eZOTIkezfv5+QkBBOnz5dqf+Plcms1KmUGgi8DxiABSIyo8j2McC7QGzuqrkissCCcWqV7fPBxde1uRs6ToCMFFh6X/HtIQ9BuwhIjodvHi68beyaMt/y1KlTrFixgvnz59OhQweWLVvG1q1bWbVqFW+99RY//PBDofbz58+na9eu+Pn5MXv2bHbs2EHdunWZO3cus2bNIiwsrMz33LNnD1u3bqVGjRrMnz+f2rVrs3v3btLT0+natSv9+/cv9ujv+PHjfPbZZ3Tt2pVx48bx8ccfM2XKFCZMmMCmTZto3rw5I0aMAKB+/fosWLCAWbNmsXr16jLjqcrKPHNQShmAj4A7gdbAg0qp1iaafi0iIbkv6yaGE78Yf2C1Ks3Pz4+goCDs7Oxo06YNffv2RSlFUFAQ0dHRxdo3aNCA1157jd69ezN79mzq1q1b7vcMDw+nRo0aAPzyyy98+eWXhISE0KlTJ+Lj4zl58mSxfXx8fOjatSsAI0eOZOvWrRw7dgw/Pz/8/f1RSjFy5Mhyx1LVmXPm0BE4JSJRAEqp5cBQ4EhlBlZhV87CsvugyxQY8Ka1o7Edpf2ld3Qpfburh1lnCkU5OTnlf21nZ5e/bGdnR1ZWlsl9Dh48iIeHR6k3IO3t7cnJyQEo9pzf1dU1/2sR4cMPP2TAgAGlxln0kWDecnV/FGzOPQcv4FyB5ZjcdUUNV0odUEqtVEr5mDqQUupRpVSkUioyLi6uAuGa4fQm47/tRlXO8TWr2bVrF+vWrWPfvn3MmjWLM2fOAFCrVi2uXbuW387X15c9e/YA8O2335Z4vAEDBvDf//6XzMxMAE6cOEFycnKxdmfPnmX79u0AfPXVV3Tr1o1WrVpx5syZ/HsKX331lWU+ZBViTnIwlR6L1rP/CfAVkWBgI/CFqQOJyHwRCRORsHr1yqw1UTGnN0GtxlCvZeUcX7OK9PR0JkyYwMKFC2ncuDGzZ89m3LhxiAhjxozhsccey78hOX36dJ588km6d++OwWAo8Zjjx4+ndevWhIaGEhgYyMSJE02esQQEBPDFF18QHBxMQkICkyZNwtnZmfnz5zN48GC6detG06ZNK/PjW4eIlPoCugA/F1h+EXixlPYG4GpZx23fvr2Yq+kLq0tdzpedJfK2j8h0N5GPu4rk5Jj9HreaI0eOWDsEm3DmzBlp06aNtcOoMFPfZyBSyvj9FBGzzhx2A/5KKT+llCPwALCqYAOlVKMCi+HA0RtNWhVyfj+kXQXvDvD3QUhJsEoYmlYdlHlDUkSylFJTgJ8xnhUsFJHDSqnXMGagVcATSqlwIAtIAMZUYswl8wqFSdvh4hGI2Q0Jp403yzStgnx9fTl06JC1w7AKs/o5iMhaYG2Rdf8u8PWLGC83rEspaNAaDLnj1+NPg09H68akaTbKJnpImiXtKvwwGS4cBPemoOwgIcraUWmazao+yeHMFti/BNISwd4RAu+F2t7WjkrTbJZNDLwyy+lN4FjTeDMSYPin1o1H02xc9TlzOL0JfLsbzxry5GSDFO2SoWmaOapHckiIgstnoFmBYdz7l8GbDSH5kvXi0qqVgwcP0rBhw1vm6UX1SA6J5403IQsmBxdPyM4wPs7UNAt46623+N///sdbb71l7VBuiupxz8G3Gzz5R+F1Hs2M/8afhiadb35MWrWTN35i2bJlVo7k5rD95JA7+g67IidB7k1AGfTjTE2rINu/rIjZBbP8IWZP4fUGB6jTVF9WVHG33357mW0++OADAgICiIiIqPR4XnnlFWbNmlVs/cGDB2natCn//e9/y3W89evX07JlS5o3b15mLczs7GzatWuXX4vy3Llz9O7dm4CAANq0acP7779fqL2vry9BQUGEhISYVeimvGzyzMHLvQa+04z1AyIMG3nT4RLUali8Ydg4cHa/ydFp5fG///2vzDYff/wx69atM2tylvxBQ7lnkkWXKyooKIjly5fzzDPPMGnSJLP2yc7OZvLkyWzYsAFvb286dOhAeHg4rVubqpUE77//PgEBASQmJgLGuhSzZ88mNDSUa9eu0b59e/r161do/99++w1PT88b+mwlsckzh23T+uTXk3yzf25RUVcTQ8Bvnwqhuq5DVVazZk2io6MJCAhgwoQJtGnThv79+5OaaqwT+thjjxEVFUV4eDhz5sxhyZIldOzYkZCQECZOnEh2dnb+/o8//jihoaFs2bKl0PK5c8ZyJKb2BXjzzTdp2bIld9xxB8ePHy8x1vr163P48GGzP9uuXbto3rw5t912G46OjjzwwAP8+OOPJtvGxMSwZs0axo8fn7+uUaNGhIaGAsaaFQEBAcTGxprcvzLYZHIoJCWeRKlRuH9DHhG49jdkpd/8uLRyOXnyJJMnT+bw4cO4u7vnF2n55JNPaNy4Mb/99hsDBw7k66+/Ztu2bezfvx+DwcDSpUsBY53Hhx9+mH379tG0adNiy0ePHjW57549e1i+fDn79u3ju+++Y/fu3SXGOG3aNNLT0/nzzz/z13Xv3p2QkJBir40bNxIbG4uPz/W6R97e3iX+cj/11FO88847JZ7hREdHs2/fPjp16pS/TilF//79ad++PfPnzzf/P9tMNnlZUUhKPJelFm6mtp36FZYOh7HroWmXmx2ZzZi5aybHEo5Z9Jit6rbihY4vmN3ez8+PkJAQANq3b2+yhuSvv/7Knj176NDB2As2NTWV+vXr06NHD5o2bUrnztefShVdLmnfhIQE7rnnHlxcXABKLIO/fv16kpOTGTx4MIcPH84v7rJly5YSP9OKFSuKrTNVWm716tXUr1+f9u3b8/vvvxfbnpSUxPDhw3nvvfdwc7v+k75t2zYaN27MxYsX6devH61ataJHjx4lxlNetp8cfDry0750ppjaVjf3GjXhtE4OVVzBepIGgyH/sqIgEWH06NG8/fbbhdZHR0cXqg0JFFsuad/33nuvzFqQaWlpPP/886xatYrPP/+cQ4cOMWjQIMB45lCwRF2eWbNm4e3tnX9JA8ZLB1Nza2zbto1Vq1axdu1a0tLSSExMZOTIkSxZsoTMzEyGDx9OREQEw4YNK7Rf3rHq16/PPffcw65duyyaHMqsBlNZrxupBGX29qxMkVfrimx4xez3ulVUlUpQrq6uxaotvfvuuzJ9+vT85aZNm0pcXJwcPnxYmjdvLn///beIiMTHx0t0dHSx/U1Vbypp3z179khQUJCkpKRIYmKiNG/eXN59991C+7700kv561asWCGjRo0y67NlZmaKn5+fREVFSXp6ugQHB8uhQ4dK3ee3336TwYMHi4hITk6OjBo1Sp588sli7ZKSkiQxMTH/6y5dusi6deuKtavsSlBVW7bpKsUAGOyNPSf148xqoXXr1rzxxhv079+f4OBg+vXrx19//XVD+4aGhjJixAhCQkIYPnw43bt3L7Tf8ePH2bBhA0899RRgfGphbvdpe3t75s6dy4ABAwgICOD++++nTZs2+dsHDRpUahXtbdu2sXjxYjZt2pR/L2PtWmNZlb///ptu3brRtm1bOnbsyODBgxk4cKBZcZlLiZUGJoWFhUlkZKRZbX2nrSF6holJVwDe8uKjlD5MfmOx6e1L7oVrF2DS1gpGWj0dPXqUgIAAa4ehVTJT32el1B4RKbNjhG3fc8hMhYwkkqVGyW06ToD04teEmqaVzraTQ4pxVqsEapXcpkXpE5Zommaabd9zyE0Ol6WU5JCVDrF7IOniTQpK06qHapEcWh+PInnHzkKbknfsJH7BAuNw7k/7wIn11ohQ02yWzSWHuA/nXl+o1Ri6TGGvuz+xTz+dnyCSd+wk9umncQ4Mgto+YOeQPzqz0P6appXI5u45XProI1J27Sq0blBUPPb163N2/Hjs69cn6+JFnJo149JHH3EJaOzjRMb65Vz6IoqU3bupN9VklylN0wqo0skhJTOF1VGrsXNMMrldqUzypvI0uLkZE8P589g3boyhQDfTzExXHBxSbkbImlZtVOnLisycTGbumolDncL3E5ou/pKmi7+kyYjGNGnzP17o/jiekycjqal4Pj7J+O/kyfntXLoPwNHDmaaLv7TSJ9FuptJqRJhTP8KUkuo8mKtXr17k9etZsWIFAQEB9O7du8LHuxmq9JlDbafa9GnSh3UZm8nIzsDR4Ijn5MnXG6TEg4snwSdOEfv0m3jNmYNr5064dOxE7NNP5y/TYQK0MfZLL7S/lq/rjE3EXik+nqGivNxrsG1an7IbVgJTNSKys7MxGAxm1Y+obJ999hkff/xxlU8OVX5sxZaYLRK4KFB+if6l+MZFQ0QW9Jd/DPuHJG3fUWhT0vYdcunTT816j1tR0T73ZY1fKS9zjrd48WLp0KGDtG3bVh599FHJysoSEeN4i+eff15CQ0Olb9++snPnTunZs6f4+fnJjz/+KCIin3/+uYSHh8uAAQOkRYsW8sor18fPuLq6iohxnEKvXr3kwQcflICAgELbRERmzpwpgYGBEhwcLC+88IKIiMyfP1/CwsIkODhYhg0bJsnJySIiMn369GJjLkRERo8eLRMnTpRu3bqJv7+//PTTTyIikpKSIiNGjJCgoCC5//77pWPHjrJ792559dVXxdXVVVq0aCHPPvtsuf9fy+tGxlZU+eSQlZ0lrT+9XR7f+HjxjR/fLrLswbJ/EJMuiRz+0fivJiLWTw5HjhyRIUOGSEZGhoiITJo0Sb744gsREQFk7dq1IiJy9913S79+/SQjI0P2798vbdu2FRFjcmjYsKFcunRJUlJSpE2bNrJ7924RKZwcXFxcJCoqKv9987atXbtWunTpkv/LHx8fLyIily5d/xl56aWX5IMPPhCR0pPDgAEDJDs7W06cOCFeXl6Smpoqs2fPlrFjx4qIyB9//CEGgyE/vp49e+Z/XdluJDlU6csKAIOdgcyroWyL3UJcShz1XApUfEqJh8btyj7IpRPwzSgY9X3h8vWa1ZRUXwHA0dExfxBRUFAQTk5OODg4EBQUVKjOQ79+/fDwMM6iPmzYMLZu3VqslmLHjh1NlpfbuHEjY8eOza/jULduXQAOHTrEyy+/zJUrV0hKSmLAgLJ72N5///3Y2dnh7+/PbbfdxrFjx9i8eTNPPPEEAMHBwQQHB5fnv6dKqPLJASDzanuyPX9nddRqxgaOvb6h61Pg0Ry2l1Hpqabxh073kqw6pIT6CgAODg75NRbs7Ozyaz3Y2dmRlXV9FG7ROgym6jIUretQ8P1NtR8zZgw//PADbdu2ZdGiRSaLrxRVUhxl1Ymo6qr004o8klGPkHoh/HDqB+O1UJ7Oj4H/HWUfQCeHKqdv376sXLmSixeN35OEhIRC5dfMsWHDBhISEkhNTeWHH36ga9euZu/bv39/Fi5cSEpKSv77A1y7do1GjRqRmZmZX4KuLCtWrCAnJ4fTp08TFRVFy5Yt6dGjR/7+hw4d4sCBA+X6bFWBWclBKTVQKXVcKXVKKTWtlHb3KqVEKWXxOtlDmw8l6moUhy7ljqXPTDNOWJOZVvbOjjXBvgYk/W3psLQKupHaDHm6devGqFGj8msxlKc8+8CBAwkPDycsLIyQkJD8x5Svv/46nTp1yi+7Zo6WLVvSs2dP7rzzTj755BOcnZ2ZNGkSSUlJBAcH884779CxY8dyfbaqoMx6DkopA3AC6AfEALuBB0XkSJF2tYA1gCMwRURKLdZQ3noOB1/rQZ9v+hDeLJx/dfmXcTDVp33gwa/x/Ty75HoPed4LBp9OevbtXEXH+dvao8xFixYRGRnJ3LnW7Q4/ZswYhgwZwr333mvVOEpS2fUcOgKnRCQq98DLgaHAkSLtXgfeAZ41J+jyquVYi75N+7LuzDqe7/g8TinG00BcPPByv5Y/j0WeYj+c938BNepWRmjVgrX6JGhVlznJwQs4V2A5BuhUsIFSqh3gIyKrlVKVkhwA7m5+N2ui1rDp7CbuTL5iXOlSl23TOhRrWzRZmPVUQ7MZY8aMYcyYMdYOg0WLFlk7hEpjzj0HU7dc869FlFJ2wBzgH2UeSKlHlVKRSqnIuLg486PM1bFhRxq5NuL7k9/nD9fG1czZfmL3QuTCcr+npt2qzEkOMYBPgWVvoGBVzFpAIPC7Uioa6AysMnVTUkTmi0iYiITVq2dihqqyglV23N38brb/tZ1ziX+CnT04mZyxorgTP8PqZyA7s9zvq2m3InOSw27AXynlp5RyBB4AVuVtFJGrIuIpIr4i4gvsAMLLuiFZUcP8h2Gn7PjW2QCD/wPmPkuuWR8QSL5UGWFpWrVTZnIQkSxgCvAzcBT4RkQOK6VeU0qZnh6oEjV0bUgP7x58//d2MkMeMn/HvL4Oybqvg6aZw6x+DiKyVkRaiEgzEXkzd92/RWSViba9KuusIc/9Le4nIS2BX48uN3+nmg2M/+qOUDZl1apVZU5dv2jRolLnf6hKfH19uXTJNs5ebaL7dFG3N76dxjmKlfs+ZmCg6Vm0vdxrFHpiEeaWyErQHaEqIH7BApwDg4zD33Ml79hJ2qGDeBSYFboyhIeHlzh/ZZ5FixYRGBhocqq5kmRlZWFvb5M//jeNTXSfLspgZ+De5HR25iQRfTXaZJtt0/oQPWNw/mtfoitM2QOBw29usDYsr96mc2BQyTU6qVhdzujoaFq1asX48eMJDAwkIiKCjRs30rVrV/z9/dmVWwpw0aJFTJliLOs3dOhQvvzSWLBn3rx5REREsHLlSiIjI4mIiCAkJITU1NRCf50jIyPp1asXYCzY8uijj9K/f38efvhhsrOzee655+jQoQPBwcHMmzevwnEmJCRw9913ExwcTOfOnfO7S8fHx9O/f3/atWvHxIkTC3X/X7JkCR07diQkJISJEyeSnZ1NdnY2Y8aMITAwkKCgIObMmVPu/1uLMWfoZmW8bmiuzOxsiXu9roQsCpJ3dxUfRmvWMW5x5syVeaRlK4keOUqiR46S0+FD5UibQDnRu48caRMop8OH5m870rJVud//zJkzYjAY5MCBA5KdnS2hoaEyduxYycnJkR9++EGGDh0qIsah2ZMnTxYRkQsXLkizZs1k8+bN4u/vnz/MuugQ6Ly5NUVEdu/eLT179hQR47Dr0NBQSUlJERGRefPmyeuvvy4iImlpadK+fftCw7vLE+eUKVPya0r8+uuv+UPLp06dKq+++qqIiKxevVoAiYuLK3HIemRkpNxxxx3573/58uVy/98WVK2HbJuUdgXPrCx617yNH0//yNTQqTgZnMreb99SMDhC8H2VH2M1U1qNzory8/MjKMh49tGmTRv69u2LUqrY0Ow8DRo04LXXXqN37958//33+cOsyyM8PJwaNYwzpP3yyy8cOHCAlStXAnD16lVOnjxZbIi3OXFu3bqVb7/9FoA+ffoQHx/P1atX2bx5M9999x0AgwcPpk6dOkDJQ9bvuusuoqKimDp1KoMHD6Z///7l/oyWYpvJIbcD1H0NurDh9BI2/LmBIbcNKXu/fYuNfSN0cjBbXt3NvEsJz8cncfmr5XhOnpx/D+Joq4rNuZk3FBtKH5pd0MGDB/Hw8Cj1BqS9vT05OTkApKUVHphXcAi3iPDhhx+WWbPBnDjFxBil0oZuSylD1v/44w9+/vlnPvroI7755hsWLrRO5z2bvOdAzfpw/2I6BT5Ik1pNWHF8hXn7udbTTyvKIa/eZl5i8Jozh3pPPIHXnDmF7kHcrLqcu3btYt26dezbt49Zs2Zx5swZAGrVqsW1a9fnQ/X19WXPnj0A+X/NTRkwYAD//e9/ycw0dow7ceIEycnJFYqt4BDt33//HU9PT9zc3AqtX7duHZcvXwZKHrJ+6dIlcnJyGD58OK+//jp79+6tUDyWYJvJwbk2tA7Hzr0J97a4l70X93Li8omy96vZQD+tKIe8+T3SDh28XqwXcO3cCa85c0g7dLBQu8qUnp7OhAkTWLhwIY0bN2b27NmMGzcOEWHMmDE89thj+Tckp0+fzpNPPkn37t0xGAwlHnP8+PG0bt2a0NBQAgMDmThxYolnLGV55ZVXiIyMJDg4mGnTpvHFF18AMH36dDZv3kxoaCi//PILTZo0AUoesh4bG0uvXr0ICQlhzJgxJs8sbhpzbkxUxuuGbkjGnxY5uVEkM10up16WsMVhMn3b9LKP8ftMkeluIplpZr93dWXODUnN9t3IDUnbPHM48iMsGQY5Wbg7uzP4tsGsjlrN5bTLpe+X30uy/IO+NO1WY5vJISXeWNnJ0VgcNCIggvTsdL49WfL1JQDBI+Cf56G2900IUtNsm20mh+R4cPHIX/Sv40+nRp346thXZOaUMurSoQY4mi44qmlaYbaZHFLiwaXwM+6RASO5mHKRX//8teT90pPg55fgzOZKDtA2SBklAjXbdqPfXxtODh6FVvXw7oFPLR+WHi2lYrDBAbbPhXM7S25zi3B2diY+Pl4niGpKRIiPj8fZ2bnCx7DNTlBD5lCgGBVgLATzUKuHmLl7JocvHaaNZ5vi+9k7GR+D6r4OeHt7ExMTQ0Uqcmm2wdnZGW/vit9fs83k0Mj07EF3N7+bufvnsuToEt7uXsLzYd3XATBOHGNqJihNy2N7lxXZWcYxEpdOFdtU07Emdze/m/XR64lLKfwXMW8I946LBnYdOkbXGZtuVsSaZpNsLzmkxMOPj0PUbyY3P9TqIbJzsvnq2FeF1ucN4e4cFEBHrxoWnaNB06oj20wOUGLV6SZuTejTpA9fH/+alMyU4g2GfwYT/68SA9S06sF2k0ORpxUFjQscR2JGIitPrCy+0c72PrKmWYPt/aaYkRyC6wUT1iCML498SWbRUvQxkbBiDPUpo6u1pt3iqmVyABgbOJa/U/5m7Zm1hTekXobD3+OlbKPIp6ZZi+0lh8DhMHEzuJQ+01V3r+741/Fn0eFF5EjO9Q2uxsl06qkrlRmlptk820sONdyhUVswlN5FQynF2DZjOXXlFFtitlzfkFui3lMlVmaUmmbzbC85HF0NB03caDRhoN9AGrk2YuGhAmW2cp9y1EOfOWhaaWwvOUR+Bjs+Nqupg50Do9uMZu/Fvey/uN+40uAAni3IMTk/sKZpeWwvOaRehhrmVx2+p/k9uDu589nBz66vnLKbD7OHVUJwmlZ92F5ySEmAGnXMbu7i4MJDAQ/xe8zvHI0/WomBaVr1YnsDr1KvFKvlUJaIgAgWH17MJ398wvt93oftH/ORw/fA4ELtus7YVKhbtZd7DbZN62OJqDXN5thWcsjOgvSr5TpzAHBzdGNU61F8/MfHHEs4RqvEWHrb/QEiUGBOgdgrqUTPuJ4wCs61qWm3Gtu6rLAzwLOnoNPEcu8a0TqCWg61+OSPT8C1Hi4qHTKSKiFITasebCs5KAU165X7zAGMZw8jW4/k17O/ctw+92zhmq7roGklsa3LioQzsH8ZhI4C9ybl3j0iIILFRxbzSfwe5gCj/rOCLTnXC8d4udewYLCaZtvMSg5KqYHA+4ABWCAiM4psfwyYDGQDScCjInLEwrFC3HHY/A60HFih5FDbqTYjW4/kkz8+4XjD1iy+owM0v8PiYWpadVDmZYVSygB8BNwJtAYeVEq1LtJsmYgEiUgI8A7wH4tHCpCaYPy3ApcVeUYGjKSmQ03mteqiE4OmlcKcew4dgVMiEiUiGcByYGjBBiJScKCCK0Wrv1pKau4w63J0giqqtlNtIgIi2PDnBo4lHLNQYJpW/ZiTHLyAcwWWY3LXFaKUmqyUOo3xzOEJUwdSSj2qlIpUSkVWqOpx6mVQduDkVv59C3i4zcO4KQfe/2n0DR1H06ozc5KDqUEIxc4MROQjEWkGvAC8bOpAIjJfRMJEJKxevXrlixSMycHZ/YarObk5ujG+hi9bSWH3hd03dCxNq67M+S2LAXwKLHsD50tpvxy4+0aCKtGd78JTByxyqAcbdqN+Vhbv7Z6lJ3bRNBPMSQ67AX+llJ9SyhF4AFhVsIFSyr/A4mDgpOVCLMDODpxqWeRQzh7NefzyVQ4kHGHTOV2mXtOKKjM5iEgWMAX4GTgKfCMih5VSrymlwnObTVFKHVZK7QeeASrnYv73mbBviWWOVceXoUnJ+Dl78sHeD8jKybLMcTWtmjDr4l1E1opIC53KAxwAAB6KSURBVBFpJiJv5q77t4isyv36SRFpIyIhItJbRA5XSrT7l8CZLWW3M0cdX+xb3MkTfkOJuhrFT6d/ssxxNa2asK3u0ymXyz0is0TObvDQcvp2eJIgzyA+/uNj0rPTLXNsTasGbCc5ZGdCxrUb6gBlipIcngp9igvJF1h2dJlFj61ptsx2xlbkd4CyYHL46SmI3krHqZH08O7B/APzCW8WjkcNY9n7vPk18+j6DtqtxHbOHNKvgYOLZZODc224HA052Twb9ixpWWl8tP+j/M1582vmvfT8mtqtxHaSg0czeOkv47wVllLHF3Iy4dpf+NX244FWD/DtyW85nnDccu+haTbKdpJDHmXBqtF1fI3/Xo4G4LG2j1HLsRbv7H5Hd4zSbnm2kxxO/wbfTjAWmLWUIsmhtlNtJodMZteFXbpjlHbLs53k8PchOPiNsVScpdT2hg7jwaN5/qr7WtxHs9rNmB05m4zsDMu9l6bZGNtJDqmXQRlueERmIQYHGDwbmnTOX2VvZ8/zHZ7n3LVzLD261HLvpWk2xnaSQ958FZa85wDGitbXLhRadbvX7fT07sknf3zC38m6zqR2a7Kd5JBqwd6RBa15Gj7pXmz1Cx1fIFuyeWf3O5Z/T02zAbbTCcreGdybWv647k0h+SJkJIOja/5qn1o+TAiawNz9c9kWu42uXl2LdYoC3TFKq75sJzkMm1c5x817YnHlLNQPKLRpbOBYVket5s2db/Jd+Hcmk4Ce+EarrmznsqKyFHmcWZCjwZGXO7/MuWvn+OzQZ8W2a1p1ZjvJYcm9sP8ryx+3lOQA0KlRJwb5DeKzg58RfdV0G02rjmwiOTiRAac2QGKs5Q/u4gF3vApNupTY5LkOz+FscObNnW/qnpPaLcMmkkNtko1fWHi4NmB8NNrtKWgcUmITzxqeTA2dyo6/dvBTlC4Ko90abCI5uKvcCW8r41EmQFIcxO4ttcmIliNoV78dM3fNJC6lAmX1Nc3G2EZyIDc5VMaZA8DW/8CiwVDKJYOdsuPV218lPTudN3a8oS8vtGrPJpKDoKB+a3CtXzlvUMcXMlMgufQzAr/afkwOmcymc5v4+c+fKycWTasibCI57JZW8Ph2aFB0ik4LKeOJRUGjWo8i0COQt3e+TUKaBUeIaloVYxPJodLl9by8/GeZTe3t7Hmt62skZiQyY9eM0hvrSw/NhtlEcnjEsBYWDam8N6jT1Dji84J5s2n51/FnYvBE1p1Zh33NAlX4s7OMpfN/fgnmdoCjuXP/nNsNK8bCyQ3GNppmA2yi+3QzFQuXTlTeGzjUgHs/A5/OZbfN9UjQI2w6u4nDWd9xKXUcnjnAsvshdg/pYs/OnADmLTlKtFsttg1OgKjf4PB3ULMhDHwbAodV3ufRNAuwiTOHOiqp8p5U5GlzD7g1Mru5g50DM7rPwM6QQfeFj7FtxhDSYg7wlv1knP4ZTY/Xt7L0rWnGorRB98I/jsP9i6G2F3z7CBz4phI/jKbdOJtIDu4kQ41K6uNQ0LE1sHmW2c1vc7+Nf3Z+Dvuaxzl71wicH1nDP19+y/R8nvZO0DocRv8ETbtCeqIFA9c0y7ON5KAsP5mNSWe2wO8zIDnevPanNvLA2SN08+rG7BNLiXKrV/Y+jq7w8I/G8nRg/ntp2k1mE8nhmDQptXuzxbSLMJaqP7ii7LZntsCyB1BRv/F62Au42Lswbcs0Ls6fR/KOnYWaJu/YSfyCBddX5NXBPL8P3m8Lfyy34IfQNMuwieTwdOZk6Pl8pRw77sO51xcaBkGjEOOEvaW5eBSWRxjn0hj9E57uvrx6+6scTTjKWueTxD79dH6CCI47RezTT+McGFT8OJ4tjUlv1VS4eKzicWtaJVDW6gYcFhYmkZGRZrX1nbaG6BmDKyWOo60CcOnQIX+5plsUHvUOcv5cTzIz3Iu1NxhSaei1BVQOF2J7kJ3lkr9tbquzrPeO552NDfHddwH7+vVJ++sCri38MbiZLoxrZ0ijsc9vZGW6cCG2O+bm65Tduwk4drR8H1bTAKXUHhEJK6td1T9zuHaBLY5PwpEfb8rbpSR5k55eG4Mh0+R2B6dElF02F//qXCgxAEw44Y3vNWfe6BFHRiMPss6fJ97ZrcTEAJCT7UzCpSCcnK/g5n7aop9F026EWf0clFIDgfcBA7BARGYU2f4MMB7IAuKAcSJSdndDcyRfwscurlJ7GzZd/GWxdQ2KrshINnaUcnCG9Gs0NvVEAvjwajTT5w4nKyGeBpMmcvWzJXhOnoxr504lByAC3zxMndo+1Bn4llkxH20VUHYjTbsBZZ45KKUMwEfAnUBr4EGlVNFBDvuAMBEJBlYClivZnDe7diUN1/acPNn0howU472Fgyvh65HwTjM4ttq4rYTEAFDv6N8896Pi3aE5LO2azVsdRhW6B2GSUnDfIjAzMZQat6ZZiDlnDh2BUyISBaCUWg4MBY7kNRCR3wq03wGMtFiEqbmDmyrpUWa9qVNMb1jQFy7mfsSaDSF0FNRrVebx0g4dxO+DjwhkA4sOLyLF72G8Jswh7dDB0s8e8p5gxOyB+JPQ9oGKxa1pFmJOcvACzhVYjgFK+SnnEWCdqQ1KqUeBRwGaNGliXoR5Zw43oxNUQX1ehj//B62GgE8nsDPv9ozHeGP/heeyQzgQd4Aj2d+Q0OYBfEpLDAVtm2Mcg+HTCer6VTR6Tbth5iQHU1NMmbwBoJQaCYQBPU1tF5H5wHwwPq0wK8JajdiQHUq/m9EJqqBWg42vCnIyODG752wGfTucAcvGkhL9OIgTUMZcF3e+CyfbwW9vwvAFptto2k1gzp/DGMCnwLI3cL5oI6XUHcBLQLiIpFsmPKDFACZkPguOLmW3rWJ83HyY1/89HGrEMWzAVs68PYjoGYON4y1K4tYIujxu7Ij11x83L1hNK8Kc5LAb8FdK+SmlHIEHgFUFGyil2gHzMCaGi5YP03bd7nU7z7R/hg1/bmD+gfnm7dT1SeM9lo2vVGpsmlaaMi8rRCRLKTUF+Bnjo8yFInJYKfUaECkiq4B3gZrACmWc6PasiIRbJMIVY1nmcAyonE5QN8PDrR/mWMIx5u6fS8u6Lcvewbk29J0OaVcgJ8fs+x2aZklm9XMQkbXA2iLr/l3g6zssHNd1SX9jp2y7opJSiuldphN1NYppW6Zh5zix7J3CxlZ+YJpWiqr/Jyn1MlekprWjuGHO9s683/t9nA3O1PBZRHyqGaMxc3LgwArj0wtNu8mqfnJISeCKuJbdzgY0dG3Ih30+RNlf44lNT5CaVcqNSQDEWDZ/3fOQbbo7t6ZVlqqdHESMZw6U3CPR1gTVCyIt9gEOXjrIi1teJDsnu+TGdga44xVIiIK9xbt4a1plqtrJIScbQh5iT46/tSOxqKykNjzf4Xl+Pfsrs/fMLr2xf39jh6gt/4GsjJsToKZR1ZODwR7ueo8NOWWOLrU5I1uPZGTASBYfWczSo0tLbqgU9HgeEmPgj0qYZVzTSmAT1aerq2fDnuV80nlm7ppJHac6DLptkOmGzftC0H1Qs5Jm/NI0E6r2mUM1Z7AzMLPHTNo3aM9LW19ic8xm0w2VMnalbnnnzQ1Qu6XpMwcr8HKvge+0NfnLjevch3/rZJ75/Rnm9ZtH+wbtTe+Yfs1YbzJs3PVRnJpWSXRysIKig658p61hdb9PGL1uNFN+ncLCAQsJ8DBRzOX0Jlj7rLFrddC9Nyla7ValLyuqiLrOdfm0/6fUdKzJYxsfI+pKVPFGre4y1pTYPMvYQUrTKpFODlVIQ9eGzO83H4Vi3M/jiicIOzvo8RzEHYVjP1knSO2WoZNDFeNX24+0c48Sl5TOXd9GcNu/P8d32hq6zthkbNDmHvBoDv/3rp7FW6tUOjlUAXk3KPNehqyGrBq+hHo1nWga+CW/vdj6eg0IOwN0/we41LleJUvTKoG+IVkFlFQV6rMBnzHu53E88vMjKMeHr29o+yCEPHSTotNuVfrMoQpr5t6MhQMWki3ZuDT9hKPxuZPYqNzKfVdjjRWyNa0S6ORQxTVzb8YXA78AsWfcz+PY+/de44acHFg0CNb8w7oBatWWTg42wLe2LynRk/Cs4cnEDRPZGrvV+OSi02Pw5zZjlWxNszCdHGyEZLmzaOAi/Gr7MXXTVNZHr4fQ0eBaz9jvQdMsTCcHG+JRw4PPBnxGsGcwz//f83x5ciXSeTKc/hVi91g7PK2a0cnBxtRyrMW8fvO4o+kdvBv5LjMMiWS7eMK53dYOTatm9KNMG1F0sBb0xql+KstYyVKXdrj97sv/OlstPK0a0snBRpjuC3EXS48uZeaumVyx/4BLsZ54elW/wjiadejLChsXERDBnN5zcHA6z4j1ozh8ZKW1Q9KqCZ0cqoG+TfqSFT0ee2XH6N2vsTpqtbVD0qoBfVlRTSSn+/FV89H848invLjlRU4knGDlxmDOX7k+bWmpE/hqWhE6OVQjdTtPYf6uBcx0duXzw5+TVas5e55cgEcND4AiNzQ1rXT6sqI6cXTBofuzvHw+hlfbPoHBJZr7f7r/epdrTSsHnRyqm/aj4akDDAuZQEr04zjbOzPu53F8fuhzQNd/0Mynk0N1Y+8ELnVBhLrprnw95Gv6NOnDf/b8hxreX5CQlmDtCDUboe85VBNFO0l96PABK2uco6bhPmb3nM2yY8t4e8e7DPtxGG90e4NuXt2sGK1mC3RyqCaKPYU4bgdfPQC7PkV1eZyIgAj+tTyVOqFrmLRxEhEBETwV+hTO9s7WCVir8sy6rFBKDVRKHVdKnVJKTTOxvYdSaq9SKksppWumVwUtBkKzvvD7DEi+BEBOeiOWD1nOyICRLD26lAfXPMjh+MNWDlSrqspMDkopA/ARcCfQGnhQKdW6SLOzwBhgmaUD1CpIKRj4NmQmw6bXAeOlR8uXNvLf7wJJOTuO0/EXiVgTwQd7PyAjW0/SqxVmzmVFR+CUiEQBKKWWA0OBI3kNRCQ6d5ueTKEqqdcSOj4KJ9ZDRkqRS4/B+P7Tm5GD/uDTg5+y6ewmXu/6OkH1gqwWrla1mHNZ4QWcK7Ack7uu3JRSjyqlIpVSkXFxcRU5hFZevV+CSdvB0aX4thwX3uj2Bh/3/ZikzCRGrhvJu7vfJSUz5ebHqVU55iQHZWJdhR6Yi8h8EQkTkbB69epV5BBaeTnVBAdnyEiBvwvfX8h7wjFqbiKn9k0iLSGML498SfgP4fz656+InhfjlmbOZUUM4FNg2Rs4XznhaJXm20fg/D54fLtxrk1MDQMfzm2vfoxbu4089ftT9PTuyYudXsSrZoVOFDUbZ86Zw27AXynlp5RyBB4AVlVuWJrF9XwekuNg7XOlNstJbcrXQ77m2bBn2XVhF0N/GMrcfXP1pcYtqMzkICJZwBTgZ+Ao8I2IHFZKvaaUCgdQSnVQSsUA9wHzlFL6+VhV07gd9HwBDq6AQ9+W2tTBzoHRbUaz6u5V9PHpw7wD87jrh7tYHbVaX2rcQszq5yAia0WkhYg0E5E3c9f9W0RW5X69W0S8RcRVRDxEpE1lBq1VULdnwCsMVj8DiWVfGTZ0bcg7Pd/hi4Ff4FnDkxe3vMjIdSPZd3HfTQhWszY9tuJWYrCHe+ZBo7aQnWn2bqENQvlq8Fe8dvtrnE86z8PrHmbqpqmcunyqEoPVrE13n77VeDaH0ebfMuo6Y9P1SXxxpHGd53nizhgWHlrI8J+GE94snMfbPk6jmo0qJ17NavSZw60qOR5WjoMLh0ptFnsllegZg/Nf5y/nMCF4AuuGrWNkwEjWRK1h0PeDeH376/yV9NdNCl67GXRyuFXlZBqn0Vt2v1n3H4pyd3bnuQ7PsXbYWob7D+e7U98x6PtBvLHjDZ0kqgmdHG5VtRrCQ99A2lVjgki/VqHDNHRtyMudX2btPWsZ1nwY3578lkHfDeLlrS8TdSXKwkFrN5O+53AraxQM9y2CZSNgxVh4cHmxuhBe7jXMO1TNRvyry78YHzSeRYcX8d3J7/jx9I/09unNI0GP0LZe20r6EFpl0cnhVuffDwbPNg7tToypUHXqwjctwcu9Ez8/NZGvjn3FsqPL+O3cb2SnNCEjoStZ1wIBg66EbQN0ctAgbCwEDgPn2pCTYxzurUwNqTEt76ZlHt9pa6jrXJfJIZMZ22Ysbf/zBv7++znn8hUNXBrwYKsHeX15rcr4JJoF6XsOmpFzbeO/v7wEP0yCLMvUd3BxcCHzcld+uvsnPuzzIb5uvry39z1qNn+bf275J/sv7te9LqsofeagXScCzu6w42NIjIURS64njVzFJ/Q1776Ewc5AL59e9PLpxcnLJxmy6F02ndvET1E/4V/Hn/ta3Mcgv0HUdqpd5rG0m0MnB+06paDXC+DuA6umwsKBELECanvnN7HEfQL/Ov6k/z2UTVPeY92ZdXx9/Gve2vkW7+5+l75N+nJ387vp3KgzBjvDDb+XVnE6OWjFhTwEbo3h61GwaAhMiTR2vbYwFwcXhrcYzvAWwzkaf5QfTv3AmjNrWB+9ngYuDRjkN4jBtw2mRZ0WqHLcA9EsQycHzbTbesEjG+DSCWNiEDFeahQ4i7CkAI8AAjwC+EfYP/j93O+sOr2KxUcW8/nhz2lWuxmDbxvMQN+B+Lj5lH0wzSJ0ctBKVr+V8QXwx1ew+mnjsO9Oj5kuO2cBjgZH+vv2p79vfy6nXeaX6F9Yc2YNH+z7gA/2fUBA3QDj9qb9aeLWpFJi0Ix0ctDMc1svaH4H/Poq/O8DCB0NHSeYPJOoaEeqouo412FEqxGMaDWC2KRYNv65kV+if+H9ve/z/t73aVmnJb2b9Ka3T28C6gboSw8LU9Z6jBQWFiaRkZFmtfWdtqbQc3TNiv7cbnyacWw1NAyCiZuN67MywN7R7MMU7ThlStGOUnn7KPsr2LsdxLXOUcQpmhzJoYFLA3r59KKHdw86NOxADfuKJaRbgVJqj4iEldVOnzlo5dO0i/F15Swk5VYQT0uE/7QGn47GHpe39QIP/1JvYprz1KPrjE3FzkCKdrba+0oXNsds5rezv7Hq9Cq+Pv41TgYnwhqG0d2rO10bd6WpW1N9VlEBOjloFePexPgCyEqH0Ifh1AZYnzshmsEJhs2DNvfAtQsQs9t4CeLmDa6eoBTxCxbgHBiEa+dO+YdN3rGTtEMH8Rg/vswE4uVeg9BXtgMOQH9QvWlQ7zzDul1ja+xWZuyaAUAj10Z0btSZLo270KlRJ+o6162E/5DqR19WVBFxH86l3tQp1g7jxl2OhrM74MJBaDcS6gfA4e9hxZjrbQxO4FqP5Jb/JPaNj/CaNh5XwyGSz6YSu3AHXlOG4NouEFoONHbCSoqD1ARwqmV8ObiCnenOvQV/Vs4lnmP7X9vZfn47Oy/s5FqGceRpc/fmdGjYgQ4NO9C+QfsqkSxu5vff3MsKnRyqiKOtAnDp0MHaYVQKpTJxcEjGYJ+KvUMKBvs0DIZ0riS0IiMhi4yoE9RpnsTVKGe8br+MawNj1+3YP+8gK8sVN/eT1PHIn2ANEZAce2LP9SUn2xnXWn/i4voXOTkO/H1VSE03kJFpz+mYhjgYHGjv50wOORx3zeGPukkcrJPEEfdk0g3GCdq8k51ofaUmra+40vpKTRqlOqJMTtdSeVJ27ybg2NGb8l76noNWZYg4kJHhDhnuxbYZ3MBQrxEJR89j37gR8ckdSIjOxM4ui6ws403FlOSGZGXVwE5lYWeXhbIzbpccBwCUysZgn4aD3TX8XDKwM2QhAh4OgeyISsC9znFq1Y6mKdAr05608w4kRjvx9bnm/NkonUtNk9neIIFfvOIBqJ1hoOVVV1pddaXlVVf8E11wyb71emvq5FCFNF38pbVDuOmSd+wk9umn8Xx8Epe/Wo7H408UugdRGreSNmRnodKu0tTVg57T1jBiSgO4cACS41m2YRdjOtbALSebfzywFIBN/+pJr9T9RDk4sNfZiQNOjhzwEHbVSwSMU775OtahTS1f2tQLok3jLrRsGIqLg+X6ehxtFWCxY1mKTg5VhOfkydYO4abLSwxec+bg2rkTLh07FVquMIM9uHpcX/YOM76AV9atYcz9hS9Rn8ycwsGpgTRPPE/zxPPcnxgLBgeudhzP4UuHObD+KY4kx7IrJY7V8fvg2JcooKmbLy3rtqRV0lVaujTCv2E7GjQIQdX2KXd386r4/dfJoYqoFjcjyynt0MFCicC1cye85swh7dDBG0sO5eTm7oHvf/JK2tXFy92LbdP6UBu43et2bh+3E1Li4cqfxF08zJGL+zkqqRx3cubQpUP8nBRr3PXs99TKzsE/MxN/d3+atR1Jc/fm3HZmB3Xr+qM8moF7U5P9Qari918nB81qPMaPL7bOtXOnm5oYoHifi6JD0lHK+PjV1ZN6Xu3pycP0LLA5Me0yp2J3cuJCJCcTjnEyKYY1GX+TtPOt/Da1s7Pxy8zCLzMLX/ua+LW4i6btH8HHuT4Of/4PPJoZHw0bHCrxk5aPTg5atWaJrtymenMW7L3p5lyH0GYDCW02MH+7iHAx5SKnr5wm6uIBTsUdYHP0EY7XuEqqfRacWwXnVmGHolFmJk2yMmmSmY2PQy28XRviHTIG71ZDcRWM/USskDh0ctCqNUvUnyhaBg+K994snQ/gg5f7veya1ofEjESir0Yz7NMfeXqAG39ePMC5a+dYlxZPomRAzl+w923Y+zbu9q40Tk7AKyuHxvYuNHL2pFEtbxq1HUXDxmG42zmhlB04ON/w5yxKJwdNK8Kcs40bSTpujm4E1wsmK/Eck9sXTjpX068SkxRDzDXj6/zlk8ReOsKp1EtszkoiXeIgMQ62GOcrdVb2NMhIpQH2NLB3pb6zB/VrNqJ+i7vwdPelvmMdPGs2xNHeqdxx6uSgaUXcrKrYppLQtml9qO1UmzYexeeiFhES0hK4kHyBC8l/8VfyBf66eICLcUf4Oz2BPVkpXEw7S1b6Odi+q9C+tQU8lSOeDjXNjk8nB02zkjJvhBahlMKjhgceNTxo45mbPFoXbpMjOVy+dp64zGvEpV7iYvT/cenCfuJSLxGfmUhcxkWz49PJQdOqsLKGthcd1m6n7PBw88YDaAXg3b3YPksfNa9ruE4OmlZFlFTZu7RxRabONopPMlSxCYTMSg5KqYHA+4ABWCAiM4psdwK+BNoD8cAIEYkudzSadguryC+wOQnF/KcqhZWZHJRSBuAjoB8QA+xWSq0SkSMFmj0CXBaR5kqpB4CZwIgKRaRpmtnMSSimEog5zDlz6AicEpEoAKXUcmAoUDA5DAVeyf16JTBXKaVET2WkaVZXNIGomebtZ05y8ALOFViOAYr2b81vIyJZSqmrgAdwqVBQSj0KPJq7mKSUOm5emHiqmYWPVQ14gv5MNqA6fqaW5jQyJzmYurVZ9IzAnDaIyHxgvhnvWfjgSkWaU5zClujPZBuq62cyp505E+nGYOz/mccbOF9SG6WUPVAbSDAnAE3TqiZzksNuwF8p5aeUcgQeAFYVabMKGJ379b3AJn2/QdNsW5mXFbn3EKYAP2N8lLlQRA4rpV4DIkVkFfAZsFgpdQrjGcMDFo6z3JciNkB/Jttwy34mqxWY1TStajPnskLTtFuQTg6applUpZODUmqgUuq4UuqUUmqateO5UUopH6XUb0qpo0qpw0qpJ60dk6UopQxKqX1KqdXWjsVSlFLuSqmVSqljud+zLtaO6UYppZ7O/dk7pJT6SilVYpWYKpscCnTbvhPjwNQHlVKtS9+ryssC/iEiAUBnYHI1+Ex5ngRuzqwsN8/7wHoRaQW0xcY/n1LKC3gCCBORQIwPGEp8eFBlkwMFum2LSAaQ123bZonIXyKyN/fraxh/2LysG9WNU0p5A4OBBdaOxVKUUm5AD4xP4hCRDBG5Yt2oLMIeqJHbH8mF4n2W8lXl5GCq27bN/yLlUUr5Au2AndaNxCLeA54HcqwdiAXdBsQBn+deLi1QSrlaO6gbISKxwCzgLPAXcFVEfimpfVVODmZ1ybZFSqmawLfAUyKSaO14boRSaghwUUT2WDsWC7MHQoH/ikg7IBmw6fteSqk6GM++/YDGgKtSamRJ7atycjCn27bNUUo5YEwMS0XkO2vHYwFdgXClVDTGS78+Sqkl1g3JImKAGBHJO7NbiTFZ2LI7gDMiEicimcB3wO0lNa7KycGcbts2RSmlMF7DHhWR/1g7HksQkRdFxFtEfDF+jzaJSIl/jWyFiFwAziml8kYw9qVwmQJbdBborJRyyf1Z7EspN1mrbJm4krptWzmsG9UVGAUcVErtz133TxFZa8WYtJJNBZbm/nGKAsZaOZ4bIiI7lVIrgb0Yn5zto5Su1Lr7tKZpJlXlywpN06xIJwdN00zSyUHTNJN0ctA0zSSdHDRNM0knB03TTNLJQdM0k/4fG7s+/MTFgTMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ex = np.linspace(0, xmax)\n",
    "pl.hist(\n",
    "    x, ex, density=True, fill=False, histtype=\"step\", label=\"empirical pdf\"\n",
    ")\n",
    "\n",
    "pl.plot(ex, mixture_pdf(ex[:, None]), \"--\", label=\"mixture pdf\")\n",
    "\n",
    "pl.plot(\n",
    "    ex,\n",
    "    hat_lambda * np.exp(-hat_lambda * ex),\n",
    "    label=f\"inferred $\\hat\\lambda$={hat_lambda:.3f}\",\n",
    ")\n",
    "\n",
    "pl.errorbar(\n",
    "    mu_.squeeze(-1),\n",
    "    np.exp(log_prior_) / (2 * np.pi * s2_.squeeze(-1)) ** 0.5,\n",
    "    None,\n",
    "    s2_.squeeze(-1) ** 0.5,\n",
    "    \"x\",\n",
    "    capsize=2,\n",
    "    label=\"mixture modes\",\n",
    ")\n",
    "\n",
    "pl.legend(loc=\"upper right\")\n",
    "pl.xlim(xmax=xmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Environment (conda_mxnet_p36)",
   "language": "python",
   "name": "conda_mxnet_p36"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
